L’objectif de ce TP est de vous permettre de vous familiariser avec les packages R utilisés pour la manipulation et la visualisation de variations génomiques à partir d’un fichier VCF. Le Jeux de données choisi pour ce TP est le même que celui utilisé pour le TP variant du niveau 1.

Pour ce TP nous allons avoir besoin des packages R ci-dessous :

library(vcfR)        # extraire les informations stockées dans un VCF
library(tidyverse)   # collection de packages R pour la manipulation et la visualisation des données,
library(reshape2)    # transform les data frame
library(ggplot2)     # visualisation
library(UpSetR)      # upset plot -> graphique d'intersection
library(venn)        # diagramme de venn

1 Description du jeux de donnée :

voir lien

Samples QTL*
SRR1262731 QTL-
SRR1205992 QTL+
SRR1205973 QTL+
  • QTL+ : diminution de la production en lait et une augmentation des concentrations en protéine et lipide

Pour le TP nous disposons des résultats du variant calling des 3 échantillons (en Multi-VCF annoté).

# le chemain vers le vcf
vcf_file <- "data/data1_mutation/pool_GATK_annot.vcf"

#samples <- c("SRR1262731", "SRR1205992", "SRR1205973")

2 Lire le fichier vcf

La fonction read.vcfR() permet de lire un fichier vcf/multi-vcf et retourne un objet de la class vcfR.

# pour chercher l'aide d'une fonction utiliser la fonction help() ou ?
# help(read.vcfR) 
# ?read.vcfR()
vcf <- vcfR::read.vcfR(file = vcf_file)
## Scanning file to determine attributes.
## File attributes:
##   meta lines: 43
##   header_line: 44
##   variant count: 2826
##   column count: 12
## 
Meta line 43 read in.
## All meta lines processed.
## gt matrix initialized.
## Character matrix gt created.
##   Character matrix gt rows: 2826
##   Character matrix gt cols: 12
##   skip: 0
##   nrows: 2826
##   row_num: 0
## 
Processed variant 1000
Processed variant 2000
Processed variant: 2826
## All variants processed
# pour vérifier le type ou la class d'un objet on utilise summary(), class(), is()
summary(vcf)
## Length  Class   Mode 
##      1   vcfR     S4

Les différentes sections d’une class S4 s’appellent des “slots”. La fonction slotNames() permet de lister les noms des slots de notre objet

#help(slotNames)
slotNames(vcf)
## [1] "meta" "fix"  "gt"

Pour accéder au contenu d’une section (slot) –>

@meta

## @meta est un vecteur de caractères
# afficher les 6 1er éléments du vecteur
head(vcf@meta)
## [1] "##fileformat=VCFv4.2"                                                                                                                 
## [2] "##ALT=<ID=NON_REF,Description=\"Represents any possible alternative allele not already represented at this location by REF and ALT\">"
## [3] "##FILTER=<ID=LowQual,Description=\"Low quality\">"                                                                                    
## [4] "##FORMAT=<ID=AD,Number=R,Type=Integer,Description=\"Allelic depths for the ref and alt alleles in the order listed\">"                
## [5] "##FORMAT=<ID=DP,Number=1,Type=Integer,Description=\"Approximate read depth (reads with MQ=255 or with bad mates are filtered)\">"     
## [6] "##FORMAT=<ID=GQ,Number=1,Type=Integer,Description=\"Genotype Quality\">"

@fix

## @fix est une matrix R
head(vcf@fix) %>% data.frame()

Nombre de lignes correspond au nombre de variants == 2826 variants.

@gt

## @gt est une matrix R
head(vcf@gt) %>% data.frame()

Chaque ligne correspond aux informations de génotypage par variant et par échantillon.

2.1 Quelques fonctions utiles du package vcfR

2.1.1 extract.indels()

Séparer les SNVs des INDELs

# extraire les SNVs
vcf.snv   <- extract.indels(vcf, return.indels = FALSE)

Nombre de SNV == 2355.

# extraire les InDels
vcf.indel <- extract.indels(vcf, return.indels = TRUE)

Nombre de InDels == 471.

2.1.2 tidy vcfR

Convertir l’objet vcfR en format compatible avec les fonctions tidyverse.

2.1.2.1 vcf_field_names()

Cette fonction extrait les metadata sous format tabulé.

# extraire les champs INFO du metadata
stringr::str_subset(vcf.snv@meta, pattern = "INFO=")[1:3]
## [1] "##INFO=<ID=AC,Number=A,Type=Integer,Description=\"Allele count in genotypes, for each ALT allele, in the same order as listed\">"
## [2] "##INFO=<ID=AF,Number=A,Type=Float,Description=\"Allele Frequency, for each ALT allele, in the same order as listed\">"           
## [3] "##INFO=<ID=AN,Number=1,Type=Integer,Description=\"Total number of alleles in called genotypes\">"
## convertir la section metadata @meta en data frame type tibble
#help(vcf_field_names)
vcf_field_names(vcf.snv, tag = "INFO")[1:3,]

La même chose pour le champs FORMAT :

# extraire les champs INFO du metadata
stringr::str_subset(vcf.snv@meta, pattern = "FORMAT=")[3:5]
## [1] "##FORMAT=<ID=GQ,Number=1,Type=Integer,Description=\"Genotype Quality\">"                             
## [2] "##FORMAT=<ID=GT,Number=1,Type=String,Description=\"Genotype\">"                                      
## [3] "##FORMAT=<ID=MIN_DP,Number=1,Type=Integer,Description=\"Minimum DP observed within the GVCF block\">"
## convertir la section metadata @meta en data frame type tibble
#help(vcf_field_names)
vcf_field_names(vcf.snv, tag = "FORMAT")[3:5,]

2.1.2.2 extract_info_tidy()

Cette fonction permet d’extraire la colonne INFO du vcf sous format tabulé. Cela facilite la manipulation des ces info par la suite.

# extraire les champs INFO de la section @fix
vcf.snv@fix[,"INFO"][1]
## [1] "AC=4;AF=0.667;AN=6;BaseQRankSum=2.97;DP=40;ExcessHet=3.9794;FS=0.000;MLEAC=4;MLEAF=0.667;MQ=59.16;MQRankSum=0.00;QD=20.23;ReadPosRankSum=1.65;SOR=0.534;ANN=A|5_prime_UTR_variant|MODIFIER|ENSBTAE00000385309|GENE_ENSBTAE00000385309|transcript|transcript:ENSBTAT00000055044|protein_coding|1/16|c.-59T>A|||||72287|,A|intragenic_variant|MODIFIER|ABCG2|ENSBTAG00000017704|gene_variant|ENSBTAG00000017704|||n.37913396T>A||||||,A|non_coding_transcript_variant|MODIFIER|gene:ENSBTAG00000017704|gene:ENSBTAG00000017704|transcript|ENSBTAT00000055044|protein_coding||||||||WARNING_TRANSCRIPT_NO_START_CODON"
## convertir le champ INFO en data frame type tibble
#  on peut choisir les 
#extract_info_tidy(vcf.snv, info_fields = c("AC", "AN", "MQ"))
extract_info_tidy(vcf.snv)

2.1.2.3 extract_gt_tidy()

Cette fonction permet d’extraire chaque metriques et scores de la colonne FORMAT sous format tabulé.

# 6 1ere lignes de la section @gt
vcf.snv@gt[1:6,] %>% data.frame()
## convertir la section @gt en data frame type tibble
#  on peut choisir les champs à garder
#extract_gt_tidy(vcf.snv, format_fields = c("AC", "AN", "MQ"))
extract_gt_tidy(vcf.snv)

2.1.2.4 vcfR2tidy()

Cette fonction convertit toutes les info stockées dans l’objet vcfR en data frame de type tibble.

# la fonction "vcfR2tidy" prend en entrée un objet vcfR
# pour connaitre le type d'un objet on fait la commande : is(vcf.snv)
# on peut specifier les champs que on veux garder (ex ci-dessous) 
# vcf.tidy.list <- vcfR2tidy(vcf.snv, format_fields = c("GT", "AD", "DP"), info_fields = c("AC", "AN", "MQ")) 
vcf.tidy.list <- vcfR2tidy(vcf.snv)

Cette fonction retourne une liste de 3 tibbles fix, gt, meta

names(vcf.tidy.list)
## [1] "fix"  "gt"   "meta"
head(vcf.tidy.list$meta)
head(vcf.tidy.list$fix)
head(vcf.tidy.list$gt)

Pour encore facilter la manipulation de nos données :

On va fusionner les deux tibble $fix et $gt :

# "ChromKey", "POS" sont les colonnes communes
vcf.tidy <- dplyr::full_join(vcf.tidy.list$fix , vcf.tidy.list$gt, by = c("ChromKey", "POS"))
#head(vcf.tidy)

Le champs ID renseigne un identifiant unique pour chaque variant. S’il n’est pas renseigné il est recommandé d’en créer un pour faciliter la manipulation des tableaux.

## exmple :
# CHROM <- "Chr1"
# POS <- 10000
# REF <- "A"
# ALT <- "T"
# ID  <- paste(CHROM, POS, REF, ALT, sep="_") # résultats : Chr1_10000_A_T

# Pour chaque variant si ID n'est pas null, créer le ID (comme dans l'ex.), sinon garder le ID.
vcf.tidy <- dplyr::mutate(vcf.tidy, ID = dplyr::if_else(is.na(ID), paste(CHROM, POS, REF, ALT, sep="_"), ID))
head(vcf.tidy)

l’information lié au nombre de reads par allele est stocké dans la colonne gt_AD (Allele depth) sous forme de : (ref AD,alt AD) ou (ref AD,alt1 AD,alt2 AD) dans les cas multi-alléliques.

Dans le cas des variants multi-alléliques les valeurs de la colonne ALT contiennent une “,” pour séparer les différents allèles alternatifs.

# les variants multi-alléliques
#stringr::str_detect(vcf.tidy, string = ALT, pattern = ",")
#dplyr::filter(vcf.tidy, stringr::str_detect(string = ALT, pattern = ",")) %>% select(ID, ALT, gt_AD) %>% head
vcf.tidy.multiAlt <- dplyr::filter(vcf.tidy, str_detect(string = ALT, pattern = ","))
vcf.tidy.multiAlt %>% select(ID, ALT, gt_AD) %>% head

Pour le reste du TP on va s’intéresser uniquement aux variations avec un allele alternatif

# garder les variations  avec un allele alternatif
filter(vcf.tidy, !str_detect(string = ALT, pattern = ",")) %>% select(ID, ALT, gt_AD) %>% head
vcf.tidy <- dplyr::filter(vcf.tidy, !str_detect(string = ALT, pattern = ","))

On va séparer la colonne gt_AD en deux colonnes ref AD et alt AD. Pour cela on utilise la fonction “separate” du package “tidyr”.

# separarer la colonne gt_AD en deux colonne ref AD et alt AD
vcf.tidy <- tidyr::separate(data = vcf.tidy, col = gt_AD, c("gt_ref.AD", "gt_alt.AD"), sep = ",", remove=FALSE)
select(vcf.tidy, ID, Indiv, gt_AD, gt_ref.AD, gt_alt.AD) %>% head()

On va ajouter une colonne avec l’allèle ratio de chaque variant.

#vcf.tidy <- dplyr::mutate(vcf.tidy, gt_AR = gt_alt.AD/gt_DP)
#vcf.tidy <- dplyr::mutate(vcf.tidy, gt_AR = if_else(gt_DP == 0, 0, gt_alt.AD/gt_DP))

vcf.tidy <- dplyr::mutate(vcf.tidy, gt_AR = if_else(gt_DP == 0, 0, as.numeric(gt_alt.AD)/gt_DP))

select(vcf.tidy, ID, Indiv, gt_AD, gt_ref.AD, gt_alt.AD, gt_AR) %>% head()

3 Contrôle de qualité et nettoyage des donnnées :

p <- ggplot(data=vcf.tidy) + 
  geom_density(aes(QUAL, color=Indiv)) + 
  xlab("score de qualité") + ylab("densité") + xlim(c(0,100))
p

p <- ggplot(data=vcf.tidy) + 
  geom_density(aes(gt_DP, color=Indiv)) + 
  xlab("Variant position coverage") + ylab("densité") + xlim(c(0,100))
p

p <- ggplot(data=vcf.tidy) + 
  geom_density(aes(as.numeric(gt_alt.AD), color=Indiv)) +
  xlab("Alternative allele depth") + ylab("densité") + xlim(c(0,50))
p

p <- ggplot(data=vcf.tidy) + 
  geom_density(aes(gt_AR, color=Indiv)) + 
  xlab("Alternative allele ratio") + ylab("densité")
p

La visualisation des différents métriques nous permet de voir la qualité de nos données mais aussi de définir les filtres et seuils à appliquer.

Ci-dessous on vous donne quelques exemples de filtres mais il existe plein d’autres filtres qualité

# filtrer les variations avec un score qualité inférieur à 30
vcf.tidy.flt.q <- dplyr::filter(vcf.tidy, QUAL >= 30)
# les positions couvertes à moins de 4 reads seront concidéré comme non couverte
# l'allèle alternatif doit être supporté par au moins 2 reads.
vcf.tidy.flt <- dplyr::mutate(vcf.tidy.flt.q, gt_DP = if_else(gt_DP < 4, 0, as.numeric(gt_DP)),
                                            gt_alt.AD = if_else(gt_alt.AD < 2, 0, as.numeric(gt_alt.AD)))

# on peut aussi appliquer des filtre sur l'allèle ratio

Petit exercice :

Dans le cours Filtre et annotation (Filtre et annotation) de EBAIIn1 nous avons appliqué un filtre qualité avec la fonction VariantFiltration et SelectVariants de la suite GATK. Appliquer les mêmes filtres sur l’objet vcf.tidy.flt.q en utilisant la fonction “filter” de “dplyr”.

#

4 Statistiques discriptives

4.1 Croisement des variants entre les 3 échantillons

On besoin d’extraire une matrice des allèles ratio (pour quoi AR?) avec les variants en lignes et les échantillons en colonnes

# transformer la matrice de données en format "wide" -> ID en lignes et Indiv en colonne
vcf.mat.AR <- reshape2::dcast(vcf.tidy.flt, ID~Indiv, value.var="gt_AR")
head(vcf.mat.AR)
# retirer la colonne ID -> pas utile
vcf.mat.AR <- dplyr::select(vcf.mat.AR, -ID)

# remplacer les NA par des 0 
vcf.mat.AR <- dplyr::mutate_at(vcf.mat.AR, .vars = names(vcf.mat.AR),
                                           .funs = function(x){dplyr::if_else(is.na(x), 0, as.double(x))})
head(vcf.mat.AR)

On peut comparer les variants des 3 échantiilons avec un diagramme de venn.

# la fonction venn necessite une marice binaire : absence, présence.
vcf.mat.bin <- dplyr::mutate_at(vcf.mat.AR, .vars = names(vcf.mat.AR), 
                                            .funs = function(x){dplyr::if_else(x > 0, 1, 0)})
head(vcf.mat.bin)
# dessiner le diagramme de venn
venn::venn(vcf.mat.bin, zcolor = "style")

ou un upset plot :

UpSetR::upset(vcf.mat.bin, sets=names(vcf.mat.AR),
              order.by = "freq", nintersects = NA,
              mainbar.y.label = "SNV Intersections", 
              sets.x.label = "Number of SNV")

4.2 Clustering hiérarchique des accessions.

# calcul la distance entre les lignes de la matrice
dist  <- dist(t(vcf.mat.AR), method = "euclidean")
# clustering hiérarchique
hc    <- hclust(dist)
# arranger la structure du dendrogramme
dend  <- as.dendrogram(hc)

# dessiner le dendrogramme
plot(dend, main = "Sample clustering using all variants")

5 Recherche de la mutation responsable

La mutation recherchée :

  • se situe dans un gène codant pour une protéine,
  • variation missence (changement de l’acide aminé),
  • présente seulement sur le locus QTL- (SRR1262731)

Appliquer les filtres nécessaires pour cerner la mutation d’interêt.

Pour garder que les variants sur des gènes codant et les variants annotés missences.

#L'information de l'annotation est decrite dans la colonne "ANN"
# vcf.tidy.flt.q$ANN[1]

# 1) garder que les variants situées dans des gènes codants et annotées comme "missense_variant"
vcf.tidy.flt.ann <- dplyr::filter(vcf.tidy.flt, grepl("protein_coding"    ,ANN), 
                                                grepl("missense_variant"  ,ANN))
vcf.tidy.flt.ann %>% head()

Pour sélectionner que les variants qui portent l’allèle référent sur les échantillons QTL+ et l’aalèle alternatif sur l’echantillons QTL-.

# extraire une matrice avec variants (position) en lignes et échantillons en colonne.
# avec l'allele ratio  (value.var)
vcf.mat.AR.ann <-  reshape2::dcast(vcf.tidy.flt.ann, ID~Indiv, value.var="gt_AR")
vcf.mat.AR.ann
# filtre
mutations <- filter(vcf.mat.AR.ann, SRR1262731 !=  0 & SRR1205973 == 0 & SRR1205992 == 0)
mutations
resultats <- dplyr::filter(vcf.tidy, ID == mutations$ID)

resultats